% *******************************************************************************************
% m_eval_1945_50_robustness.m
% ********************************************************************************************
clear all; 
close all;
format shortG;
% ************************ MAIN ROBUSTNESS ********************************
% 1. Low value of discount factor (rho/sigma=8; rho=0.5)
% 2. Low theta after 1960 (theta=0.3) 
% 3. High cost share for floor space (gammaH=0.2) 
% *************************************************************************
for c=1:4
% ==================== LOAD FUNDAMENTAL AND DATA  =========================
load parameters.mat;                     % basic parameters 
load calib_productivity_amenity.mat;     % inverted fundamentals every 5 years from 1955
load calib_setup.mat;                    % Info data every 5 years from 1945
load commuting_cost.mat;  
load calib_auxiliary.mat;                % Info used in this calibration 

Rdata45 = Rdata_all(:,1); Rdata50 = Rdata_all(:,2); 
Ldata45 = Ldata_all(:,1); Ldata50 = Ldata_all(:,2); 
N = size(Ldata_all,1); 
% =================== PARAMETERS ==========================================
alpha=alphaest;   beta=betaest; 
% =========== COMPUTE LOWER BOUNDS for Theta 1950 & SET Theta 1950 ======== 
Rthetalb = (Rdata50 - Rdata45 < 0);   Rthetalb = (1 - Rdata50./Rdata45).*Rthetalb;  Rthetalb = max(Rthetalb); 
Lthetalb = (Ldata50 - Ldata45 < 0);   Lthetalb = (1 - Ldata50./Ldata45).*Lthetalb;  Lthetalb = max(Lthetalb);
thetalb = ceil(100.*max([Rthetalb; Lthetalb]))./100;

theta50 = 0.90;
if min(theta50 - thetalb)<0
    disp('>>> CHECK THETA VALUES (LOWER BOUND) ! <<< ');   pause
end 
theta55 = thetavec(1); 
% ===================== COUNTERFACTUAL SET ================================ 
if c==1; 
    rho=0.5; rho50=0.5; sigma=rho/8; sigma50=rho50/8; 
end 

if c==2; 
    theta55=0.7; 
end 

if c==3; 
    gammaH=0.2; 
end 
% ==================== COMPUTE COMMUTING PATTERN 1945 =====================
options0 = optimset('Display','none','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-5,'TolX',1e-5);
K45=Kappa_mat(:,:,1).^(-rho/sigma); 
Y0=ones(N,1);     
comsyst = @(y)fcomeval2(y,K45,Rdata45,Ldata45,rho,sigma);
[y_fsolve, yfval_fsolve]=fsolve(comsyst, Y0, options0); Y45=abs(y_fsolve); Y45=Y45./geomean(Y45);  
ytemp=K45.*(repmat(Y45,1,N).^(rho/sigma)); 
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
Lcom45=ytemp.*repmat(Rdata45',N,1); 
%  ================ INVERSION OF (X,Y) FOR 1950 =============== 
options0 = optimset('Display','none','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-6,'TolX',1e-6);
X0=ones(N,1);  
Y0=ones(N,1);  
Rsyst = @(x)fReval2(x,K50,Q50,Rdata45,Rdata50,Ldata45,Ldata50,rho50,sigma50,theta50,mu,gammaL,gammaH);
Lsyst = @(y)fLeval2(y,K50,Q50,Rdata45,Rdata50,Ldata45,Ldata50,rho50,sigma50,theta50,mu,gammaL,gammaH);  
[x_fsolve, xfval_fsolve]=fsolve(Rsyst, X0, options0);   
[y_fsolve, yfval_fsolve]=fsolve(Lsyst, Y0, options0); 
X50=abs(x_fsolve); X50=X50./geomean(X50);
Y50=abs(y_fsolve); Y50=Y50./geomean(Y50);
% ************************************************
% *********** FORWARD LOOKING MODEL **************
% ************************************************
% =========== DECOMPOSE (X,Y) INTO FUNDAMENTALS AND AGGLOMERATION =========
% compute fundamental productivity and amenities in 1950
b50 = X50.*((Rdata50./Area).^(-beta)).*(X55.^(-rho*(1-theta55))); 
a50 = Y50.*((Ldata50./Area).^(-alpha)).*(Y55.^(-rho*(1-theta55))); 
% ========== COMPUTE STRUCTURAL RESIDUALS (FORWARD LOOKING) ================
l_ave_b = sum(log(bb),2)./size(bb,2); % average amenities over time (1955-75) 
l_ave_a = sum(log(aa),2)./size(aa,2); % average productivity over time (1955-75)
l_bb50 = log(b50) - l_ave_b;         % (log) structural error in amenities in 1950 (Forward Looking)
l_aa50 = log(a50) - l_ave_a;         % (log) structural error in productivity in 1950 (Forward Looking)

% Note: Since we have included the time-fixed fundamental for 1950, this
% is actually fundamentals + 1950 year fixed effects. 
l_b50_ave = sum(log(b50))./N; 
l_a50_ave = sum(log(a50))./N; 
bb_err50 = exp(log(b50)-l_b50_ave); 
aa_err50 = exp(log(a50)-l_a50_ave); 
XX50 = exp(l_ave_b).*((Rdata50./Area).^beta).*(X55.^(rho*(1-theta55)));   
YY50 = exp(l_ave_a).*((Ldata50./Area).^alpha).*(Y55.^(rho*(1-theta55))); 

% ********************************************************************
% ******** SIMULATE FORWARD LOOKING MODEL ****************************
% ******** WITHOUT STRUCTURAL RESIDUALS ******************************
% ** NOTE: HERE WE USE PROBABILITIES THAT IS BASED ON XX50, YY50 *****
% ********************************************************************
XX_temp = K50.*((repmat(Q50',N,1).^(-mu)).*repmat(XX50',N,1)).^(rho50/sigma50); 
XX_temp = XX_temp./repmat(sum(XX_temp,2),1,N); 
YY_temp = K50.*((repmat(Q50,1,N).^(-gammaH/gammaL)).*(repmat(YY50,1,N).^(1/gammaL))).^(rho50/sigma50); 
YY_temp = YY_temp./repmat(sum(YY_temp,1),N,1);

R50_p=Rdata50; 
L50_p=Ldata50; 
maxit=1e+6; 
ii=0; 
update=0.05;
while ii < maxit
    Ldif_p = L50_p-(1-theta50).*Ldata45;
    Rdif_p = R50_p-(1-theta50).*Rdata45; 
    
    R50_n = (1-theta50).*Rdata45+sum(XX_temp.*repmat(Ldif_p,1,N),1)'; 
    L50_n = (1-theta50).*Ldata45+sum(YY_temp.*repmat(Rdif_p',N,1),2);
    R50_p_r = round(R50_p*1e+6)./1e+6; R50_n_r = round(R50_n*1e+6)./1e+6; 
    L50_p_r = round(L50_p*1e+6)./1e+6; L50_n_r = round(L50_n*1e+6)./1e+6; 

    if min(R50_p_r == R50_n_r)==1 && min(L50_p_r == L50_n_r)==1
        break; 
    else
        ii=ii+1; 
        R50_p = (1-update).*R50_p+update.*R50_n; 
        L50_p = (1-update).*L50_p+update.*L50_n;
        if rem(ii,200)==1
        disp(['Largest Error is:',num2str(ii), ' ',num2str(max(abs(R50_p_r-R50_n_r))), ' ',num2str(max(abs(L50_p_r-L50_n_r)))]);
        end
    end
end 

% ==================  SAVE RESULTS ========================================
if c==1
    R50_smallrho = R50_p; L50_smallrho = L50_p; 
end 

if c==2
    R50_largetheta = R50_p; L50_largetheta = L50_p; 
end 

if c==3
    R50_largegammaH = R50_p; L50_largegammaH = L50_p; 
end 

end 

result=array2table([ID, Rdata45, Rdata50, Ldata45,Ldata50, R50_p, L50_p, R50_smallrho, L50_smallrho,...
    R50_largetheta, L50_largetheta, R50_largegammaH, L50_largegammaH, cbddist, Area],'VariableNames',{'geocode',...
    'pop45','pop50','emp45','emp50','pop50_base','emp50_base','pop50_smallrho','emp50_smallrho','pop50_largetheta','emp50_largetheta',...
    'pop50_largegammaH','emp50_largegammaH','cbddist', 'area'}); 
writetable(result,'../../output/quant/output_evaluations_1945_50_robustness.csv');

% ************************************************************************
% ************************************************************************
% **************************  FUNCTIONS **********************************
% ************************************************************************
% ************************************************************************
function[Reval] = fReval2(X,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH)   % eq D.16. X corresponds to Theta

N = size(DRpre,1); 
X=abs(X);  % avoid negative values
X=X./geomean(X);
xtemp=K.*((repmat(Q',N,1).^(-mu)).*repmat(X',N,1)).^(rho/sigma);
xtemp=xtemp./repmat(sum(xtemp,2),1,N); 
Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 
Reval = Rdif - sum(xtemp.*repmat(Ldif,1,N),1)'; 
Reval = abs(Reval); 

end

function[Leval] = fLeval2(Y,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH) % eq D.17. Y corersponds to Omega

N = size(DRpre,1); 
Y=abs(Y);  % avoid negative values
Y=Y./geomean(Y);
ytemp= K.*((repmat(Q,1,N).^(-gammaH/gammaL)).*(repmat(Y,1,N).^(1/gammaL))).^(rho/sigma);
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 
Leval = Ldif - sum(ytemp.*repmat(Rdif',N,1),2); 
Leval = abs(Leval); 

end

function[comeval] = fcomeval2(Y,K,DR,DL,rho,sigma) % eq D.17. Y corersponds to Omega

N = size(DR,1); 
Y=abs(Y);  % avoid negative values
Y=Y./geomean(Y);
ytemp=K.*(repmat(Y,1,N).^(rho/sigma)); 
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
comeval = DL - sum(ytemp.*repmat(DR',N,1),2); 
comeval = abs(comeval); 

end
